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We study the surface barrier for magnetic field penetration in mesoscopic samples of both type I and 
type II superconductors. Our results are obtained from numerical simulations of the time-dependent 
Ginzburg-Landau equations. We calculate the dependence of the first field for flux penetration 
(Hp) with the Ginzburg-Landau parameter (k) observing an increase of Hp with decreasing k for 
a superconductor- insulator boundary condition ((V — iA)'ii\„ = 0) while for a superconductor- 
normal boundary condition (approximated by the limiting case of ^'|s =0) Hp has a smaller value 
independent of k and proportional to He. We study the magnetization curves and penetration fields 
at different sample sizes and for square and thin film geometries. For small mesoscopic samples we 
study the peaks and discontinuous jumps found in the magnetization as a function of magnetic field. 
To interpret these jumps we consider that vortices located inside the sample induce a reinforcement 
of the surface barrier at fields greater than the first penetration field Hpi. This leads to multiple 
penetration fields Hpi = Hpi, Hp2, Hp3, . . . for vortex entrance in mesoscopic samples. We study the 
dependence with sample size of the penetration fields Hpi. We explain these multiple penetration 
fields extending the usual Bean-Livingston analysis by considering the effect of vortices inside the 
superconductor and the finite size of the sample. 
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I. INTRODUCTION 

In the last years there has been an important exper- 
imental and theoretical interest, in the study of vortex 
physics in a mesoscopic scale. Em The smallness of these 
systems imply that the sample geometry and the inter- 
action between the vortices and the sample surface be- 
come important. The interaction between vortices and 
the surface manifests itself fundamentally in the exis- 
tence of a|-surface barrier, first studied by Bean and 
Livingston,El which delays the vortex penetration and 
generates metastable states. If the surface effects are ig- 
nored, the penetration of magnetic field is energetically 
favorable at the first critical field Hd- However the en- 
ergy barrier of the surface prevents the vortex entry until 
a higher field Hp at which the barrier vanishes. Hp, also 
known as the superheating field, is associated with the 
peak in the magnetization curves and is strongly influ- 
enced by the presence of surface irregularities. 

The surface barrier has attracted a renewed interest re- 
cently in the study of mesoscopic superconductors. For 
example, Enomoto and OkadaQ by means of numeri- 
cal simulations of the time-dependent Ginzburg-Landau 
equations (TDGL), studied the influence of temperature 
and surface irregularities on the surface barrier. Sonin 
and TraitoQ showed that the presence of the surface bar- 
rier affects the entry and exit of vortices influencing the 
surface resistance. They found a surface-induced sup- 
pression of the ac losses. 

One important line of research of mesoacntiippSjaper- 
conductors are the superconducting disks .El'D'OOilJ The 
study of small superconducting disks was started by Buis- 
son et ala for disks with radius ~ 7 fim. Recent ad- 
vances in the microfabrication technology and measure- 



ment techniques now allow the fabrication and study of 
superconducting disks with sizes comparable ta the co- 
herence length ^, with radius as small as 0.3 /xm.cl Most of 
the studies were done in Al disks, a material with k ~ 0.3, 
however for small samples the effective penetration depth 
A = X^/d increases for decreasing disk thickness (d) re- 
sulting in effective k values in the type II region that can 
be studied theoretically using the equilibrium Ginzburg- 
Landau equations. In this regime the Al disk can develop 
Abrikosov multivortex stateso and depending on the ra- 
dius (i?) and thickness of the disk it is cossible to observe 
first or second order phase transitions,!! by increasing the 
disk sizes the second order reversible phase transition ob- 
served for small disk radius is replaced by a first order 
transition. There is also an intermediate regime where 
jumps in the magnetization appear associated with the 
vortex entrance. Other interesting phenomena have been 
studied for mesoscopic Al disks, for example in Ref. |l^ 
hysteresis in the magnetization curves was observed ex- 
perimentally and explained in terms of the presence of 
a "Bean-Livingston" surface barrier and in Ref. [ij the 
behavior of the third critical field Hc3 was investigated 
for different sample sizes and geometries. 

The time dependent Gipaburg-Landau (TDGL) equa- 
tions have been proposed^ as a time dependent gener- 
alization of the mean field approach ot^the Ginzburg- 
Landau theory. Gorkov and Eliasberglla obtained the 
TDGL equations from the microscopic BCS theory in the 
gapless case. In the last years, numerical simulations of 
the time dependent Ginzburg-Landau (TDGL) equations 
have been suscesfully used to study the magiietic proper- 
ties-and fiux dynamics in superconductors. Ej E3 Frah et 
al.tli and Liu et al&B sinmlated the TDGL equations fpt 
K = 0.3 - 20, Kato et alH and Machida and KaburakiB 
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for K = 2, Aranson et studied vortex dynamics in the 
K = oo limit , and Vicente- AlvareZ|-e|i al. studied the dy- 
namics of d-wave superconductorsB3 Only few studies of 
the TDGL equations been done in superconductors 
with type I behavior ,E-3't3 possibly because the theory is 
better to describe a superconductor near a second order 
phase transition at a temperature near Te- 
la this paper, we present a numerical simulation of the 
TDGL equations to study the surface barrier in meso- 
scopic samples for k = 0.15 — 2. We neglect demagneti- 
zation effects and therefore we assume that the sample is 
infinite in the direction of the external magnetic field (the 
z direction). We consider square samples that are meso- 
scopic in the xy plane (perpendicular to the magnetic 
field) with linear sizes of 5 — 30 A, with A the penetration 
depth. 

The paper is organized as follows. In Sec. II we present 
the TDGL equations with their discretized form in fi- 
nite differences, and we discuss their possible boundary 
conditions. In the Sec. Ill of the paper we study the de- 
pendence of the penetration field Hp with the Ginzburg- 
Landau parameter k {k — A/^), exploring both the type 
I (k < I/V2) and the type II region (k > f/\/2) for large 
samples. We study the effects of the surface barrier from 
a comparison of two types of boundary conditions, (i) 
The superconductor-insulator (S-I) boundary condition: 
consisting in the vanishing of the superconducting cur- 
rent perpendicular to the boundary (Jg -6 = 0). In this 
case we find an increase of Hp with decreasing k. (ii) 
The superconductor-normal (S-N) boundary condition: 
approximated as the vanishing of the superconducting 
order parameter at the boundary (\I'|5 = 0). A different 
behavior is observed for the S-N boundary condition, the 
field Hp is independent of k and nearly equal to He- In 
the Sec. IV we study the surface barrier in mesoscopic su- 
perconductors. In particular, in Sec.IVA we study mag- 
netization curves in type II superconductors at differ- 
ent sample dimensions in the region where the transition 
from a macroscopic to a mesoscopic behavior takes place. 
In Sec.IVB we show that the discontinuities that appear 
in the magnetization curves of mesoscopic samples can 
be explained by considering that the vortices that are 
inside the sample induce a reinforcement of the surface 
barrier at fields greater than the first penetration field. 
In this way, it is possible to define a second, third, fourth, 
etc. penetration fields which are a consequence of the in- 
teraction between vortices and the surface currents. We 
study the sample size dependence of the first, second and 
third penetration fields and we show that for sufficiently 
large sample sizes the known macroscopic behavior is re- 
covered, i.e. a continuous magnetization curve appears 
since Hp^ Hp2 Hp. Finally in Sec.V we give a 
summary of our results and conclusions. 



II. MODEL AND DYNAMICS 



A. TDGL equations 

Our numerical simulations are carried out using 
the time-dependent Ginzburg-Landau equations comple- 
mented with the appropriate MaxweU equations. In the 
zero-electric potential gauge we have:t§t3 

^ = i[(V - ^A)H + (I - r)(l - (1) 
= (1 -r)Im[**(V- iA)^'] - K^v X V X A (2) 
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where ^I^ and A are the order parameter and ?.i«ctor poten- 
tial respectively and T is the temperaturcEJ Equations 
(1) and (2) are in their dimensionless form. Lengths 
have been scaled in units of ^(0), times in units of 
io = 47rcrAVc2 ^ i{{)f/riD, A in units of i?c2(0)C(0) 
and temperatures in units of Tc. rj is proportional to the 
ratio of characteristic times for 5* and A, r/ ~ tij,/to ~ 

/{AnaK^D), with = ^"^ /D, where a is the quasiparti- 
cle conductivity and D is the electron diffusion constant. 
For superconductors with magnetic impurities we have 
D — / {ASiTK^a), and therefore 77 = 12 in this case. 

We have used the standard finite differenpe discretiza- 
tion scheme to solve equations (1) and (2)B3 The order 
parameter and vector potentials are defined at the nodes 
of a rectangular mesh {f — (/, J)), and the link variables 
U i^j = exp^—iKhp^Apj j) (fj. — x,y) are introduced in 



order to maintain the gauge invariance under discretiza- 
tion. 

In our simulations we have assumed a sample that has 
a square/rectangular shape in the x, y direction with di- 
mensions Lx X Ly and it is infinite in the z direction. We 
apply the magnetic field parallel to the z direction, the 
symmetry of the problem then implies for all mesh points 
Ai^j = {A,j,i,j,Ayi^j,0) and B/^j = (0, O.B^/^j), where 
S /,j = (V X A), = {d^A^i^j - dyA^i^j). 

In this geometry the discretized form of equations (1) 
and (2) are: 



dt T^K (Aa;)2 



-I- 



T . 



+ (i-^)(i-i*/,jr)*/,./ 

dA^i^j ^ r Ini[t/,/,j^l^j^/+i,j] 
dt ■ tJ Ax 

- ( -f. ) 

Ay 

dt ^ tJ Ay 
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(3) 



(4) 



(5) 



where Aa; and Ay are the mesh widths of the spatial 
discretization. 
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B. Boundary conditions 

The dynamical equations must be completemented 
with the appropriate boundary conditions for both the 
order parameter and the vector potential. 

The boundary conditions for the vector potentials 
■^fj.i,J 3're obtained by making 

B = V X Ha 

at the sample surface (where Ha is the applied magnetic 
field). 

The boundary conditions for the order parameter de- 
pend sensibly on th£-physical nature of the boundary. In 
general is given hyxB 

(M)|„ = (V-zI)|„vI/ = i (6) 



where & is a surface extrapolation length which embodies 
the surface suppression (or enhancement if 6 < 0) of the 
superconducting order parameter. 

For the boundary between a superconductor and an 
insulator (or the vacuuin)-|We have b ~ ^^(0)/a, with a 
the interatomic distanccEj For low temperature super- 
conductors b is huge (~ 1cm), and the superconductor- 
insulator (S-I) boundary condition is usually approxi- 
mated with the limit 6 — s- cx): 

(V-a)|„* = (7) 

This boundary condition implies that the perpendicu- 
lar component of the superconducting current is equal to 
zero at the surface ( Js|„ — 0). This is the most frequently 
used boundary condition because it also minimizes the 
free energy at the sample surface. More precisely, this 
boundary condition is valid for superconductors with in- 
terfaces for which b ^ £,{T). 

For the boundary between a superconductor and a nor- 
mal metal we have b ^ -jff- 4r ^1°'' , with TV the local den- 
sity of states at the Fermi surface, Nn the bulk density 
of states, Tj the transmission coefficient at the bounds, 
ary and ^at is the coherence length in the normal state.E£l 
Typically {Nn /N)Tj ^ 1 and b is small when compared 
to £,{T). Therefore, the superconducting- normal (S-N) 
boundary condition can be approximated by 6 « 0, giv- 
ing: 

*|s = 0. (8) 

The case of 6 — > is— also found for a ferromagnet- 
superconductor surfaceEj Moreover, the condition \l/|s = 
is similar to having a high density of defects at the in- 
terface and therefore a highly defective surface is also 
represented by (8). It is also interesting to note that 
for high Tc superconductors even the superconductor- 
insulator boundary is better approximated by (8) since 
b ~ ^(0)^/a ^ £,{T) in a large range of temperatures, 



due to the smallness of ^(0) in this case. The bound- 
ary condition of (8) completely suppresses the currents 

at the boundary {J^ — jt — Q) and maximizes the sur- 
face Helmholtz free energy, which becomes equal to the 
free energy of a normal metal. 

In the study of the surface barrier in Sec. Ill we will 
compare these two conceptually different boundary con- 
ditions of the order parameter. They represent the two 
limiting cases of Eq. (6) and we will call them, in short, 
the "S-I" boundary condition (Eq. (7)) and the "S-N" 
boundary condition (Eq. (8)). A previous discussion of 
these two itypes of boundary conditions was done by Buis- 
son et aLlH, where the equilibrium solutions and eigen- 
values of the linearized Ginzburg-Landau equations were 
compared to study the behavior near Hc2- However, 
for the vortex nucleation process at low magnetic fields 
the nonlinear terms of the Ginzburg-Landau equations 
should be considered. 



III. SURFACE BARRIER IN MACROSCOPIC 
SAMPLES 

When the magnetic field H is increased starting from 
= in a finite sample, the Meissner state is destroyed 
at a magnetic field Hp which is typically higher than 
Hci- This is due to the presence of a surface barrier 
for vortex entrance in finite samples. The surface bar- 
rier in macroscopic samples is sometimes known as the 
"Bean-Livingston barrier" |Since it was first studied by 
Bean and Livingston (BL)eI. In the BL work the sur- 
face barrier was obtained in the London approximation 
and for the ideal case of a semi-infinite superconductor, 
therefore vortex core nucleation effects and geometrical 
effects were both neglected. With this simplification, the 
BL work was able to identify one of the most relevant 
causes of the surface barrier: the screening currents near 
the surface circulate in the opposite direction to the su- 
perconducting currents around a vortex. This is typically 
viewed as a competition between the repulsion between 
the vortex and the surface shielding current and the at- 
tractive force between a vortex inside the sample and its 
image. This argument, which is based on the London 
model, yields the Bean-Livingston result for the penc- 
tration field, i/London _ which is independent of 

k.BEj However, in addition to overcoming the barrier in- 
duced by the shielding Meissner currents, a vortex pene- 
tration event has also to deplet the superconducting order 
parameter j^*! at the surface, ijesulting in a higher pene- 
tration field Hp{K) > i^London |3a|-|35| rpj^-g ^^^^^ ^q^^^ 

be sensitive to the boundary conditions of the Ginzburg- 
Landau equations. In this section we will discuss the 
effect of boundary conditions on the surface barrier in 
macroscopic samples as a function of k. 
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FIG. 1. Magnetization curves obtained using the S-I 
boundary condition: (a) k = 0.15, (b) k = 0.8 and (c) n — 2. 
(Inserted in the figure are the size of the superconducting 
region used in the simulation, A is the penetration length). 

We start our study of the surface barrier by first an- 
alyzing the magnetization curves in large samples. The 
magnetization curves for the S-I boundary condition are 
summarized in Fig.l and for the S-N boundary condition 
are summarized in Fig. 2. In both cases the curves were 
obtained initializing the variables to a perfect Meissner 
state [^'(/, J, i = 0) = 1 and A{I, J,t = 0) = 0] and in- 
creasing the magnetic field at subsequent steps, usually 
with AH — Q.05Hc2- We take as the initial condition at 
a magnetic field H the final state of \E'(/, J) and A{I, J) 
of the previous magnetic field value H — AH. In this way 
we mimic the experimental procedure of increasing the 
magnetic field in a magnetization measurement. For each 
magnetic field we calculate the magnetization M taking 
time averages of the time dependent magnetization M{t): 

47rM(0 = ^5:i?.(/,J,<)-i?a 

^^x^V J J 

At *^ 

M^-—-J2M{t), (9) 
i Off 



FIG. 2. Magnetization curves obtained using the S-N 
boundary condition: (a) k — 0.2, (b) k = 0.7 and (c) k = 2. 

where we start the average at a time after a steady 
state was reached. In the simulations we have taken 
T = 0.5, 77 = 12 and we used a mesh of 120 x 120 
points. In order to make efficient calculations we have 
chosen the time step (At) and the spatial discretization 
(Ax and Ay) depending on the value of n. For example, 
for K = 0.15 we used Ax ^ Ay = 0.05 and At = 0.0025; 
for K = 0.8, Aa; = Ay = 0.3 and At = 0.0025 and for 
K = 2, Ax = Ay = 0.5 and At = 0.015. Since the 
TDGL eqyuit|i^iip,|«jtfisidered here represent a mean field 
dynamics ,0 lid E!lEill23 the effect of themial noise fluctua- 
tions beyond mean field are neglected, Efl which is correct 
for low Tc superconductors. 

Figure 1(a) shows the case of a type I superconduc- 
tor with K = 0.15. We can see that the TDGL equa- 
tions reproduce the basic phenomenology of type I su- 
perconductivity characterized by a first-order magnetic 
transition. In this case the superconductivity disap- 
pears abruptly and there is no surface superconductiv- 
ity. The field profile is described by a Meissner state 
and H = Ha exp(— x/A). In this simulation the interme- 
diate state structures typicaLof type I superconductors 
are not found at equilibriumcJ since the long range in- 
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teractions hetfueea currents and demagnetization effects 
are neglected H3'l3 Results for a type II superconductor 
with K ~ 2 are shown in Figure 1(c). In this case the 
superconductor is in the Meissner state until a pene- 
tration field Hp is reached. At Hp some vortices enter 
the sample and a peak in the magnetization curve is ob- 
served. Above Hp the magnetization increases due to 
the penetration of vortices in the sample, until the es- 
tablishment of surface superconductivity for fields in the 
range HcsiT) > H > Hc2{T). The criterion for surface 
superconductivity that we use is the existence of super- 
conductivity in the surface (i.e. |^| ^ in a contour 
around the surface of width equal to the discretization 
length, Ax = 0.5^ in this case) and the exact vanishing 
of superconductivity in the bulk (i.e. I^fj = inside the 
region enclosed by the surface contour). 

The S-N boundary condition leads to a different mag- 
netic behavior as can be observed in Figure 2 where we 
have used similar parameters as those reported in Figure 
1. We observe the following differences: (i) The magne- 
tization is smaller for the same external magnetic field 
(there are less vortices). Since in this case the super- 
conducting order parameter vanishes at the surface, the 
Meissner shielding currents are nucleated at a distance 
of a few ^ inside the sample, instead of being right at 
the boundary. Therefore less vortices can stay inside the 
sample for a given magnetic field when compared with the 
S-I boundary. For example, for k = 2, we find that the 
shielding distance S between the vortices and the sample 
surface is Ssn ~ 5A for the S-N interface while for the S-I 
interface we have 6si ~ 3A. In the general case, when the 
boundary condition is described by Eq.(6) with a finite 
b, one can expect that the shielding distance will have 
an intermediate value Ssn > S{b) > 6si- (ii) There is no 
surface superconductivity above Hc2- This is an obvious 
and direct consequence of the S-N boundary condition 
that enforces ^'[5 = 0. (iii) The first penetration field 
Hp is smaller and therefore the surface barrier is lower, 
we will discuss this in detail in the following paragraphs. 

In order to understand the difference in the surface 
barrier, let us study the dynamics of the first vortex pen- 
etration just at H = Hp. Figure 3 shows the dynamic 
evolution of the order parameter near a small region close 
to an S-I interface (Figure 3 left) and close to an S-N in- 
terface (Figure 3 right) . If the interface is of the S-I type 
the order parameter at the boundary is different from 
zero in the Meissner state. When the condition for vor- 
tex entrance is fulfilled the order parameter at a bound- 
ary point has to decrease until reaching zero. Therefore 
there is an intermediate time interval when 5* = in a 
point at the boundary. Just in this moment a vortex can 
enter the sample and afterwards the order parameter at 
the boundary increases again and returns to a non-zero 
value. It is interesting to note that this process is always 
necessary for vortex penetration in S-I interfaces, there 
is always the need of an intermediate \[' = state at the 
boundary. On the contrary, a smooth entrance of vortices 
is observed for an S-N interface (Figure 3 right). 




FIG. 3. Time evolution of the spatial pattern of the order 
parameter in a small region around the boundary were a vor- 
tex entrance is taking place. For the S-I condition (left) the 
magnetic field is H — 0.199-H'c2 (0), and for the S-N condi- 
tion (right) H — 0.169-ffc2 (0). In both cases H is just above 
the first penetration field and k = 2. Gray scale from black 
(1*1 = 0) to white (1*1 = 1). 

For this boundary condition the order parameter at the 
interface is already zero. A small deformation of the re- 
gion of 4* = allows for the penetration of a vortex. Since 
there is no need of depressing the order parameter at the 
boundary, the surface barrier is much smaller in the S-N 
case. A similar dynamical behavior for vortex entrance 
would appear in an S-I interface with defects at the sur- 
face. At the defects the order parameter is depressed and 
therefore ^E* = is already established at some boundary 
points, which are preferred points for vortex entrance. 

From the magnetization curves we can obtain the first 
field for flux penetration Hp as a function of k for the 
different boundary conditions, this is shown in Figure 4. 
The K dependence of the superheating field Hp has been 
previosly calculated for the case of a saaijiifinite medium 
with the Girpburg-Landau equations.t3"EJ Matricon and 
Saint JamesEll obtained Hp^n) solving the semi-infinite 
one-dimensional GL equations, in Ref . p^-p5| the stability 
of the superheated state under small fluctuations of the 
order parameter and the veptpr potential was discussed, 
and very recently VodolazovE3 analyzed the effects of sur- 
face defects on Hp. In the case of a one dimeHsional 
semi-infinite medium, the Matricon-Saint Jamea2J solu- 
tion can be obtained solving the equations: 

„ , „ = A^* + * - 
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cPA 



(10) 



with the boundary conditions: 
dm 



H = Ha and 
A = H 



dx 

and * 



at x = 

1 sX X ^ oo 



(11) 



Solving numericaUy equations (9) with the boundary con- 
dition (10) it is possible to find a relationship among 
Ha and ^(a; — 0) — ^'o, where Hp is the maximum 
value of Ha that allows a physically meaningful ^'o, i.e. 
< < 1. The results obtained in this way are rep- 
resented by the continuous-line of Figure 4. Our simula- 
tional results, on the other hand, are a numerical solution 
of the exact problem in a two dimensional square sample 
(which has geometrical effects) and are represented by 
closed circles. The Hp values reported here are the peaks 
of the simulated dc-magnetization curves and the error 
bars correspond to the discrete field step used in the the 
magnetization curves. We see in Fig. 4 that the value of 
Hp obtained for the S-I interface is always well above Hd 
in the type II region {Hci{T) = [(In K)/y2K]ffc(T)) and 
above He in the type I region, supporting the existence of 
a "surface barrier" even in the type I case. Our numerical 
simulations show that for the S-I interface Hp increases 
with decreasing k, and Hp{K) > H^, with a behavior 
in good agreement with the result for the semi-infinite 
medium obtained by Matricon-Saint James. (However, 
for smaller mesoscopic samples the value of Hp can be 
significantly larger than the Matricon-Saint James result, 
enhancing the geometrical effect of the square sample, see 
Sec. IV). For k — > our results agree with the known be- 
havior of thcLSuper heating field in type I supeconductors. 
Hp - I/Vk.EtEI For k ^ oo we obtainjhat Hp He 



in agreement with the result of Ref. |30|^32|. Figure 4 
also shows that in the type II region the Hc2 values ob- 
tained from the S-I simulations (closed squares) are close 
to the expected values {Hc2{T) = \/2kHc{T), dotted 
line). Some differences in Hc2 appear at small k near 
the type I region, since in this region the field Hp is close 
to Hc2 and a delay of the superconducting-normal tran- 
sition could be expected. 

The behavior of Hp vs k is very different when the S-N 
boundary condition is used. In this case. Hp is indepen- 
dent of K and nearly equal to TjLondon ^ ^ggg ^j^g 
open circles of Figure 4). This result shows that in the 
case of the S-N interface, when the condition ^fjs = 
is enforced, the surface barrier is only due to the surface 
shielding currents and well described by the London -ap 



proximation value of Bean and Livingston: ^London J5|,g5| 
On the other hand, in the case of the S-I interface the 
penetration field Hp is higher due to the extra contribu- 
tion needed for the vanishing of the order parameter at 
the surface in a vortex penetration event. In a type I 
superconductor, the "surface barrier" can be interpreted 
as the barrier for penetration of the normal state from 
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FIG. 4. Critical magnetic fields obtained from the magne- 
tization curves at different k values. For the S-I boundary 
condition are plotted the superheating field (Hp) (closed cir- 
cles) and the second critical field (-ffc2) (closed squares). For 
the S-N boundary condition is plotted Hp (open circles). The 
dashed and dotted curves are the expected values of He and 
Hc2 respectively. The continuous line is the result of Hp vs k, 
obtained by Matricon and Saint- James for a semi-infinite sam- 
ple. 

the boundary. In the S-N interface of a type I super- 

Hp = He, which is 



conductor there is no barrier, 



an obvious result since the boundary condition already 
enforces the normal phase — 0) at the surface. On 
the contrary, in the S-I interface of a type I superconduc- 
tor the barrier for nucleation of the normal phase at the 
boundary can be very high. Hp 3> He, as can be seen in 
Fig.4. 

In the general case, when the boundary condition of 
a superconducting sample is described by Eq.(6) with a 
finite 6, we expect that Hp(b, k) will be in between the 
two hmit cases studied here, Hp'^{b = 0, k) < Hp{b, k) < 
Hp^{b — oo,k). Ideally, most of the superconductors 
should be closer to the behavior of Hp\b = oo,k) but 
the effect of a finite b and the presence of defects at the 
surface will give a smaller value of Hp with a lower bound 
givenby i/^^(6 = 0,K). 
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FIG. 5. Magnetization curves of square samples at different 
sizes and using the S-I boundary condition for k = 2 (a) 
30A X 30A, (b) 15A x 15A and (c) lOA x lOA. 



IV. SURFACE BARRIER IN MESOSCOPIC 
SAMPLES 

A. Finite-size effects in mesoscopic type II 
superconductors 

The magnetic behavior of mesoscopic superconductors 
is different from the behavior of bulk samples. In the 
mesoscopic scale, several maxima appear in the magne- 
tization curves which are related to the vortex entrance. 
This result is quite general and appears either in thin 
films |aLjDarallel fieldgE3 or mesoscopic superconducting 
disks.l3'E£l In this section we study the magnetic behavior 
of superconducting samples of different sizes, in partic- 
ular we cover the sample size region where a transition 
from a mesoscopic to a bulk behavior takes place. 

Figure 5 shows the dc-magnetic behavior of supercon- 
ducting square samples of different sizes with S-I bound- 
ary condition. The behavior of Fig. 5(a) is typical of the 
macroscopic samples that we have studied in the pre- 
vious section, however if we decrease the sample size to 
the mesoscopic region the continuous behavior disappears 
and some magnetization maxima followed by discontinu- 
ous jumps appear (see Figs. 5(b) and (c)). 
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FIG. 6. Magnetization curves for the S-N boundary condi- 
tion and K = 2. We show the magnetic curves at different 
sample sizes from the macroscopic to the mesoscopic-like be- 
havior: (a) 30A X 30A, (b) 15A x 15A and (c) lOA x lOA. 

We find that the discontinuities in the magnetization cor- 
respond to the penetration of new vortices into the sam- 
ple, this will be discussed in detail in Sec.IVB. 

If we change the boundary condition and we use the 
'S-N' condition we find, as it is shown in Figure 6, a dif- 
ferent mesoscopic behavior: there is an appreciable de- 
crease in the number of magnetization maxima and the 
transition between the states of the system with a differ- 
ent number of vortices seems more continuous. In other 
words, for the same sample size, fewer vortex penetra- 
tion events are needed to arrive to the normal state. The 
decrease in the number of entrance events for the S-N 
boundary in mesoscopic samples is related to the fact 
that this boundary condition allows less vortices inside 
the sample, as discussed in the previous section. The 
fact that the S-N interface needs a larger shielding dis- 
tance, SsN > Ssi, has clearly stronger consequences in 
the magnetic behavior of a mesoscopic sample. At the 
same magnetic field the mean magnetization values of 
the Meissner state are lower in the S-N case than in the 
S-I case as can be observed comparing Figures 5(c) and 
6(c). In general a mesoscopic sample with a boundary 
corresponding to < 6 < oo in Eq.(6), will have a mag- 
netization behavior in between the two limiting cases of 
Fig. 5 and Fig. 6. 
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FIG. 7. Magnetization curves of thin films with the exter- 
nal magnetic field applied parallel to the film. The curves are 

generated for decreasing thicknesses of the film and using the 
S-I boundary condition for k = 2. (a) 15A x oo, (b) lOA x oo 
and (c) 5A X 00. 
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FIG. 8. Thin film magnetization curves with the external 
field applied in the same geometry of Figure 7 but using the 
S-N condition: (a) 15A x oo, (b) lOA x oo and (c) 5A x oo. As 
it is observed for samples sizes as low as 5A x oo, there is a 
continuous transition to the normal state. 



The changes observed in the magnetization curves 
when going from the macroscopic to the mesoscopic be- 
havior are quite general and do not depend on the sample 
geometry. We also study the magnetic behavior of a thin 
film with the field parallel to its faces. In this case there 
is only one relevant dimension, the thickness of the film 
(d). We work here in the case d> X. Figure 7 shows the 
magnetic behavior of thin films with different thickness 
obtained using the S-I boundary condition. We observe 
that in this case the discontinuities in the magnetization 
curves appear at smaller film thickness and that they are 
less important than in squares samples, possibly because 
vortices in a mesoscopic square sample are more confined 
than in thin films of the same linear size. We have also 
done numerical calculations of the magnetic behavior in 
thin films using the S-N boundary condition (see Figure 
8). In particular Figure 8(c) shows that for very small 
film thickness there is a continuous change in magnetiza- 
tion without the jagged structure observed in Fig.8(a-b). 
In this case vortex lines do not penetrate the sample and 
the superconducting state disappears gradually. This be- 
havior is also present when we use the S-I boundary con- 
dition but it is necessary to explore lower sample sizes 
than the one shown in the corresponding figures. 



B. Multiple penetration fields in mesoscopic samples 

The behavior of the magnetization in mesoscopic type 
II samples is characterized by maxima followed by dis- 
continuous jumps as a function of the external magnetic 
field. In Figure 9 we show in detail the case for a sam- 
ple of size lOA X lOA and k = 2. At the same time 
we plot the total number of vortices inside the sample, 
Ny = <f{A + j^)dl. We see that each discontinuous 

jump in M(H) corresponds to an increase of AA''„ = 4 in 
the number of vortices. These jumps occur at succesive 
magnetic fields Hpi = Hpi, Hp2, Hp:^, . . ., which are indi- 
cated in the figure. In the regions of Hpi < H < 
the number of vortices Ny is constant. Therefore the only 
penetration events occur at Hpi, when one vortex enters 
in each of the four sides of the square sample. In the re- 
gion of constant vorticity, Hpi < H < ifp(i+i), one may 
think that no vortices enter the sample because they feel 
a surface barrier which is enhanced by the repulsion force 
exerted by the vortices already inside the sample. 
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FIG. 9. Magnetization curves for a mesoscopic lOA x lOA 
square sample using the S-I boundary condition for k = 2. In 
the right scale the number of vortices A'^„ is shown. 

To analyze this effect, let us extend the simple Lon- 
don approximation model of Bean and Livingston for the 
surface barrier to the case when there are extra vortices 
inside the sample and near the boundary. In the Lon- 
don approach, the Gibbs free energy (G) of a vortex line 
located a cU&tance x from the sample surface can be cal- 
culated as:Eiy 

G=-[i7.exp(--)--^i.„(-) + i7.r-i^.] 

(12) 

where Ha is the external magnetic field and Ko is the 
modified Bessel function of the second kind. In normal- 
ized units the above expression reads: 



AttG 



[Haexp{- 



1 / 2a; In K 



(13) 



where we have used the relation Hd = {lnK)/2K^Hc2- 
The first term of Equations (12) and (13) is related with 
the repulsive interaction between the vortex line and the 
external field, the second term describes the attraction 
between the vortex and the surface currents, this term 
is usuaUji. interpreted as the interaction with an image 
vortexB'E^ and the third term is the vortex self energy. 

The above expressions are valid when there are no vor- 
tices inside the sample. If there are vortices located inside 
the sample, some additional terms due to the interaction 
between the vortices are needed. In the following we will 
assume that there is only one vortex inside the sample 
located at position / and that we are analysing the Gibbs 
free energy in the same line perpendicular to the sample 
surface. If I is small, the new vortex that is trying to en- 
ter the sample now feels two additional terms. The first 
term is the repulsive force between the vortices and the 
second is the attractive interaction with the image of the 
vortex located inside the sample. The last term is needed 
in order to take into account the perturbation of the 



O 
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FIG. 10. Gibbs free energy of a vortex line as a function 
of the distance a; (A) from the sample surface. We have used 
Equation (12) for (a) H = 0.03Jfc2, (b) H = 0.1H^2, (c) 
H = 0.24//c2 and Equation (13) for (d) H = 0.24Hc2 and (e) 
H = 0.28J/c2. 

vortex already inside the sample because of its proximity 
to the surface. Both contributions are more important in 
mesoscopic superconductors as we will show below. The 
free energy that gives the correct force expression is: 



AttG 



Ha exp( 



X' 



-^K (^^ 



Ha 



X 



-) 



In K, 

^k/-^) (14) 



A 



Equations (13) and (14) are approximate expressions and 
have the limitation that they are not valid near the sam- 
ple surface, then we will only use it for a; > ^. It is useful 
to note that in a more exact treatment we should obtain 
that G ^ when x — > because the Gibbs free energy 
of a vortex located outside the sample must be zero. 

Figures 10(a), 10(b) and 10(c) were generated using 
Equation (13) for k = 2. The free energy of the vor- 
tex depends on both the applied magnetic field and the 
distance to the surface. For H < Hd the free energy 
associated with a vortex is positive at all x values, as 
it is shown in Figure 10(a), vortices are then thermody- 
namically unstable inside the superconductor and there 
is an energy cost associated with the vortex entrance. 
The thermodynamic condition for vortex penetration is 
not fulfilled until H > Hd, when the free energy of a 
vortex located well inside the superconductor becomes 
negative, see Figure 10(b). However, Figure 10(b) also 
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shows that because of the attractive interaction between 
the vortex and the surface currents a barrier to vortex 
entrance appears nears the surface. This is the "Bean- 
Livingston" surface barrier originated by the screening 
currents. Taking into account Equation (13) it is possi- 
ble to estimate Hp as the magnetic field for which the 
expression {dG/dx)x=e, becomes negative, i.e. when the 
maximum of G{x) moves inside the region a; < i; . This 
condition is fulfiUed at Hp = ^Hc {Hp w 0.24iJc2 in 
Figure 10(c)). At H = Hp some vortices penetrate the 
sample and the free energy associated with the entrance 
of a new vortex now must be caknilated using an analo- 
gous to Equation (14), the exact expression depends on 
the number of vortices and their location inside the sam- 
ple. Considering only one vortex located at x = I = 3X 
and using Equation (14) we have calculated the free en- 
ergy just after the first vortex entrance {H > 0.24ilc2), 
the results are shown in Figure 10(d). Observe that the 
free energy changes considerably from Figure 10(c) to fig- 
ure 10(d). Now there are three relevant regions: i) near 
X = 3X, there is a region where the Gibbs free energy 
increases with increasing x, the strong repulsive inter- 
action with the vortex inside the sample dominates; ii) 
there is an intermediate region where G decreases at in- 
creasing X, this means that a vortex located in this re- 
gion will be pulled inside the sample, it is possible to 
allocate vortices in this region and there is also a region 
iii) near the sample surface that repels vortex entrance, 
G increases for increasing x. The existence of regions 
ii) and iii) means a reinforcement of the surface barrier 
induced by vortex penetration and allows magnetic field 
intervals of metastable states. This energy barrier rein- 
forcement due to vortex entrance is more important in 
mesoscopic superconductors because;, in small samples, 
vortices are confined by the potential well generated by 
the sample surface and even a vortex fixed in the center 
of the sample is very close to the surface. For example 
in an lOA x lOA square sample, vortices are constrained 
to be located around the center of the sample at a; = 5A, 
because of its interaction with the surface currents. As 
a consequence vortices stay at x 5X generating a new 
surface barrier and it is necessary an important mag- 
netic field increase to allow new vortex penetration. For 
macroscopic samples the situation is quite different, there 
is a nearly continuous vortex penetration. In macroscopic 
samples the vortices that are inside the sample are not 
confined and they do not have serious restrictions in their 
movement since they can be allocated very far from the 
surface. In this small increase of the magnetic 

field is enough to accommodate new vortex lines, as can 
be observed in Figure 5(a). 

From the analysis of Figure 9 and the previous discus- 
sion we can conclude that in mesoscopic samples there are 
preferred values of magnetic field for vortex penetration, 
in this case the process of vortex entrance is discontinu- 
ous in contrast with the continuous macroscopic regini(\ 
This behavior is a consequence of the barrier to vortex 



entrance that appears after each penetration event. In 
this way we can define a second penetration field Hp2, a 
third penetration field Hp^ and so on. We observe that 
after increasing the size of the sample the vortex pene- 
tration becomes continuous {Hps — » Hp2 ^ Hp). 

The exact delimitation of the macroscopic and meso- 
scopic regimes depends on the geometry of the sample. 
We will now study in detail the behavior of a thin film 
becaiise is a simpler case with only one significant length 
scale. The size dependence of the penetration fields 
Hp, Hp2, Hp3 obtained numerically from the TDGL equa- 
tions in a thin film using the S-I boundary condition are 
summarized in Figure 11(a). For d > 15 A a continu- 
ous vortex penetration is observed, this is the region of 
a macroscopic behavior. The mesoscopic region is lo- 
cated at 2A ^ d ^ 15A in which several penetration fields 
Hp. can be distinguished. For d < 2X there is a gradual 
transition to the normal state without vortices. Figure 
11(b) shows the sample size dependence of the first, sec- 
ond and third penetration field of thin films using an S-N 
boundary condition. For the S-N boundary the scales are 
shifted to larger sizes since the shielding length is larger 
{SsN > Ssi), as we discussed before. The macroscopic 
behavior appears at d > 18A, the mesoscopic region is lo- 
cated between 5A £ d ^ 18A and for d < 5A a continuous 
transition to the normal state appears. 

It is interesting to note that Hp is size dependent in 
the mesoscopic region. This size dependence can be ex- 
plained by considering in Equation (13) the effect of two 
surfaces that are separated by a small distance d. The 
term generated by the magnetic field H cxp{—x/ X) now 
becomes H cosh[{x — d)/X]/ cosh(d/A) due to the proxim- 
ity of the other surface, the magnetic force that pulls the 
vortex inside the sample decreases when the film thick- 
ness {d) decreases. The image term also changes because 
the vortex lines that are trying to enter the sample are 
now near both surfaces and new image lines are neces- 
sary in order to satisfy the boundary condition at both 
surfaces. The application of the image method to parallel 
surfaces gives an infinite number of images, but the net 
effect is the appearance of an attractive interaction to the 
new surface. Then, in the mesoscopic region, there is also 
a decrease in the net attractive image force. If we con- 
sider only three relevant terms in the image forces, the 
normalized force / = (Att / ^oHc2)X{SG / Sx) that feels a 
vortex line that is trying to enter when the sample is in 
the Meissner state can be estimated by: 

_ „ sinh(^) 1 2^ 
"cosh(d/2A) ^ X' 

1 ,2d-2x, 1 ,2d + 2x. , , 

+ ^''^ 

where d is the thickness of the film and we have chosen 
/ positive when it repels the vortex entrance. As we 
analyzed before. Hp is usually obtained from /|x=^ = 0. 
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FIG. 11. Penetration fields obtained from the magnetic 
curves shown in Figures 11 and 12. We plot the first, second 
and third penetration fields at difi^erent film thickness for (a) 
S-1 and (b) the S-N boundary condition. Hps 7^ Hp2 7^ Hp is 
typical of a mesoscopic behavior. For large film thickness, in 
the region of the macroscopic behavior, a continuos entrance 
of vortices is recovered {Hps Hp2 Hp). Figure (c) shows 
an estimation of Hp and Hp2 using the image method and 
considering the influence of the sample size in the field pro- 
file. 



This condition leads to the following expression for Hp-. 
[^i(|)-^i(^)+i^i(^)]cosh(^) 



K2sinh(^^ 



(16) 



11(c). As it can be observed, when we decrease the thick- 
ness of the film the repulsive magnetic force decreases 
faster than the attractive term, and the Hp value in- 
creases. In Figure 11(c) we also estimate the behavior of 
the second penetration field Hp2 as a function of the sam- 
ple size. The approximate expression used was obtained 
including in the Equation (15) the extra terms related 
with the presence of a vortex line inside the sample as 
we did in Equation (14). For simplicity we have consid- 
ered one vortex located at the middle of the sample at 
X = d/2. Under this condition Hp2 becomes: 



Hp2 = iK^A).K^-^: 



2\ ' 2A 
M+2£, cosh(^) 



^)]- 



2A 'V2sinh(^t^) 



(17) 



this expression also reproduces the Hp2 sample size de- 
pendence observed in the numerical simulations. 



V. SUMMARY 

We have presented results on the study of the magneti- 
zation curves and the surface barrier for type I and type 
II superconductors. Our results show that the strength 
of the surface barrier depends on the boundary. If the 
interface is of the S-I type, the vortices that enter find a 
higher barrier for penetration since the order parameter 
at the surface has to go through an intermediate state of 
4* = 0, while for the S-N boundary condition the vortex 
entrance occurs more smoothly since the ^ = condi- 
tion is already fulfilled by the interface. In this later case 
the surface barrier is only due to the Meissner shielding 
currents and the penetration field Hp agrees well with 
the Bean-Livingston value, while in the S-I case Hp is 
higher and dependent on n. Superconductors with more 
realistic boundary conditions should lie in between these 
two cases. 

We also characterized the reinforcement of the surface 
barrier due to the presence of vortex lines inside the sam- 
ple in mesoscopic superconductors. We show that these 
new barriers allow for the existence of metastable states 
of constant vorticity as a function of magnetic field. Each 
metastable state becomes unstable at the i-th penetra- 
tion field Hpi in which one vortex enters in each side 
of the sample and the magnetization has a discontinu- 
ous jump. We study the magnetization curves at differ- 
ent sample dinicinsions and we obtain the sample size 
dependence of the first, second and third penetration 
fields. We finally show that for sufficiently large sample 
sizes the continuous macroscopic regimen is recovered, 



i.e. H, 



p3 



p2 



Hr, 



2X 



We have evaluated this behavior of Hp{d) in Figure 
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